﻿using System;
using System.Collections.Generic;
using DynaPrecision.Miscellaneous;
namespace DynaPrecision.FastFourierTransform
{
	public class StandardComplexTransform
	{
		public static int BitReverse(int input, int n)
		{
			int result = 0, _input = input;
			for (uint i = 0; i < n; i++)
			{
				result <<= 1;
				result += _input % 2;
				_input >>= 1;

			}
			return result;
		}
		public static int OrderOf2(int n)
		{
			int result = 0, t = 2;
			while (t <= n)
			{
				result += 1;
				t <<= 1;
			}
			return result;
		}
		public static Complex<T> _obtainWk_N<T>(int k, int n)
		{
			double t = (2 * Math.PI * k) / n;
			return new Complex<T>((T)(dynamic)Math.Cos(t), (T)(dynamic)Math.Sin(-t));
		}
		public static void Transform<T>(IList<Complex<T>> from, int Sign = -1)
		{
			if (Sign != -1 && Sign != 1)
			{
				throw new System.Exception("Only 1 and -1 are allowed for Sign paramter in FFT");
			}
			int orderOf2 = OrderOf2(from.Count);
			int n = (1 << orderOf2);
			if (from.Count != n)
			{
				throw new FFTDataCountNotPowerOfTwoException();
			}
			//Bit Reverse
			for (int i = 0, j = 0; i < n / 2; i++)
			{
				j = BitReverse(i, orderOf2);
				if (i >= j)
				{
					continue;
				}
				else if (j < n / 2)
				{
					Swapper.Swap(from, i, j);
					Swapper.Swap(from, n - i - 1, n - j - 1);
				}
				else
				{
					Swapper.Swap(from, i, j);
				}
			}
			int level = 2;
			int halfLevel = level >> 1;
			int segCount = n >> 1;
			Complex<T> tempE, tempO, wk, wk2;

			while (level <= n)
			{
				for (int i = 0; i < halfLevel; i++)
				{
					wk = _obtainWk_N<T>(Sign * i, level);
					wk2 = _obtainWk_N<T>(Sign * (i + halfLevel), level);
					for (int j = 0; j < segCount; j++)
					{
						tempE = from[i + j * level];
						tempO = from[i + j * level + halfLevel];

						from[i + j * level] = wk * tempO + tempE;
						from[i + j * level + halfLevel] = wk2 * tempO + tempE;

					}
				}
				halfLevel = level;
				level <<= 1;
				segCount >>= 1;
			}
		}
		public static void TransformTo<T>(IList<Complex<T>> from, IList<Complex<T>> to, int Sign = -1)
		{
			int i = 0;
			for (; i < from.Count && i < to.Count; i++)
			{
				to[i] = from[i];
			}
			for (; i < from.Count; i++)
			{
				to.Add(from[i]);
			}
			Transform(to, Sign);
		}
		public static void PackedArrayTransform<T>(IList<T> from, int Sign = -1)
		{
			if (Sign != -1 && Sign != 1)
			{
				throw new System.Exception("Only 1 and -1 are allowed for Sign paramter in FFT");
			}
			int orderOf2 = OrderOf2(from.Count / 2);
			int n = (1 << orderOf2);
			if (from.Count != n * 2)
			{
				throw new FFTDataCountNotPowerOfTwoException();
			}
			//Bit Reverse
			for (int i = 0, j = 0; i < n / 2; i++)
			{
				j = BitReverse(i, orderOf2);
				if (i >= j) continue;
				else if (j < n / 2)
				{
					Swapper.Swap(from, 2 * i, 2 * j);
					Swapper.Swap(from, 2 * i + 1, 2 * j + 1);
					Swapper.Swap(from, (n - i - 1) * 2, (n - j - 1) * 2);
					Swapper.Swap(from, (n - i - 1) * 2 + 1, (n - j - 1) * 2 + 1);
				}
				else
				{
					Swapper.Swap(from, 2 * i, 2 * j);
					Swapper.Swap(from, 2 * i + 1, 2 * j + 1);
				}
			}

			int level = 2;
			int halfLevel = level >> 1;
			int segCount = n >> 1;
			Complex<T> tempE = new Complex<T>(), tempO = new Complex<T>(), wk, wk2;
			while (level <= n)
			{
				for (int i = 0; i < halfLevel; i++)
				{
					wk = _obtainWk_N<T>(Sign * i, level);
					wk2 = _obtainWk_N<T>(Sign * (i + halfLevel), level);

					for (int j = 0; j < segCount; j++)
					{
						tempE.Real = from[i * 2 + j * level * 2];
						tempE.Imaginary = from[i * 2 + j * level * 2 + 1];
						tempO.Real = from[i * 2 + j * level * 2 + level];
						tempO.Imaginary = from[i * 2 + j * level * 2 + level + 1];

						from[i * 2 + j * level * 2] = (wk.Real * (dynamic)tempO.Real - wk.Imaginary * (dynamic)tempO.Imaginary) + tempE.Real;
						from[i * 2 + j * level * 2 + 1] = (wk.Imaginary * (dynamic)tempO.Real + wk.Real * (dynamic)tempO.Imaginary) + tempE.Imaginary;
						from[i * 2 + j * level * 2 + level] = (wk2.Real * (dynamic)tempO.Real - wk2.Imaginary * (dynamic)tempO.Imaginary) + tempE.Real;
						from[i * 2 + j * level * 2 + level + 1] = (wk2.Imaginary * (dynamic)tempO.Real + wk2.Real * (dynamic)tempO.Imaginary) + tempE.Imaginary;
					}
				}
				halfLevel = level;
				level <<= 1;
				segCount >>= 1;
			}
		}
	}
}